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Granular flows due to simultaneous vertical and horizontal excitations of a flat-bottomed cylin- 
drical pan are investigated using event-driven molecular dynamics simulations. In agreement with 
recent experimental results, we observe a transition from a solid-like state, to a fluidized state in 
which circulatory flow occurs simultaneously in the radial and tangential directions. By going be- 
yond the range of conditions explored experimentally, we find that each of these circulations reverse 
their direction as a function of the control parameters of the motion. We numerically evaluate 
the dynamical phase diagram for this system and show, using a simple model, that the solid-fluid 
transition can be understood in terms of a critical value of the radial acceleration of the pan bot- 
tom; and that the circulation reversals are controlled by the phase shift relating the horizontal and 
vertical components of the vibrations. We also discuss the crucial role played by the geometry of 
the boundary conditions, and point out a relationship of the circulation observed here and the flows 
generated in vibratory conveyors. 

PACS numbers: 47.57.Gc,45.70.-n 



I. INTRODUCTION 



In a recent experimental work, Sistla, et al. [l[ observed 
a novel granular circulation occurring in a flat-bottomed 
cylindrical container (or "pan" ) subjected to coupled hor- 
izontal and vertical vibrations. The apparatus studied 
was a commonly available industrial device used mainly 
for sieving and polishing of granular particles. Of par- 
ticular interest was the appearance of a "toroidal" cir- 
culation, in which granular flow occurred simultaneously 
in the tangential direction (i.e. a circulation in the hor- 
izontal plane about the vertical axis), as well as radially 
(i.e. particles flowed outward along the bottom surface 
and then back to the center of the pan along the top sur- 
face of the granular bed). In this mode, the trajectory of 
a single particle would therefore approximate a helix on 
the surface of a torus. Among several other phenomena, 
it was noted that the direction of the tangential flow was, 
under most of the conditions explored, opposite to that 
of the orbital motion of the center of the pan. However, 
under some conditions, the tangential flow would reverse 
direction, and move in the same direction as the pan cen- 
ter. The authors noted that the phase shift relating the 
vertical and horizontal vibrations seemed to control this 
reversal. Whatever their origin, these phenomena indi- 
cate a non-trivial interplay of the horizontal and vertical 
vibrations, and their interaction with the granular bed. 
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Since the bottom surface of the pan in the apparatus was 
a rough wire screen, the authors of Ref. [l| suggested 
that cntrainmcnt of particles to the motion of the bot- 
tom surface might play an important role in establishing 
the direction of the observed granular circulations. 

Motivated by the complexity of phenomena observed in 
Ref. , we present here a computer simulation study of 
a granular matter system enclosed in the same container 
geometry as in Ref. []], and subjected to a similar form of 
horizontal and vertical vibration. As we will show below, 
we are able to reproduce the major dynamical modes of 
the granular bed observed experimentally. We also ex- 
plore the dynamical behavior beyond the range studied 
experimentally, and discover an entire family of circu- 
lation reversal transitions, affecting both the tangential 
and radial circulations. By analyzing the interaction of 
the particle bed with the rough pan bottom, we are able 
to show how particle entrainment, along with the phase 
shift relating the horizontal and vertical vibrations, to- 
gether conspire to control the direction of the observed 
circulations. 

Although most previous work on granular circulation 
has focussed on either purely vertical orpurely horizon- 
tal vibrations (see e.g. Refs. @, S S H B ) , a growing 
body of research is emerging related to granular flows 
generated by simultaneous horizontal and vertical vibra- 
tions and/or shear [g, M, liil III! • This is significant since 
most industrial devices (e.g. seives, grinders and vibra- 
tory conveyors) generate both vertical and horizontal ex- 
citations. Understanding how the fundamental physical 
mechanisms of granular flow relate to practical applica- 
tions therefore depends on the detailed examination of 
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1: Geometry and parameters describing the pan motion. 



FIG. ±: vjreomeiry anu parame 
All symbols are denned in the 



text. 



systems, such as the one studied here, having excitations 
in more than one dimension. As in our findings, these 
earlier works on simultaneous vertical and horizontal vi- 
brations also point out the significance of the phase shift 
between the two types of excitation as being an impor- 
tant control parameter. 

The cylindical geometry of our system (with the cylin- 
der axis oriented vertically) is also a distinguishing as- 
pect of the present work. A few recent works have stud- 
ied granular flow in annular troughs, and have also ob- 
served tangential circulations similar to those studied 
here 0, U|. However, Ref. pj and the present study 
are, to our knowledge, the first works to address radial 
circulation in this geometry. 



II. 



DESCRIPTION OF THE PAN MOTION 



The shape of the flat-bottomed cylindrical pan and its 
motion are illustrated in Fig. ^ The forcing motion we 
impose on the pan in our simulations is based on the 
motion observed in the experimental apparatus used in 
Ref. p] . This motion is complex and is described in detail 
in this section. 

We assume that the pan is a rigid body. Let n be 
the normal vector to the bottom of the pan, rooted at 
the center of the bottom surface. The motion of the pan 
can then be described as a superposition of two motions: 
(i) a precession of n about the vertical axis with angular 
frequency to, and with n tilted at a constant angle 7 with 
respect to the vertical; and (ii), a circular motion of the 
center of the pan in the horizontal plane, of amplitude 
A c and also at frequency lo. Note that the direction of 
rotation of both motions is always the same, and in our 
simulations, chosen to be counter-clockwise. 

In this motion, the time evolution of the point 
(x c , y c , z c ) at the center of the pan bottom has only hor- 
izontal components which are non-zero, and thus can be 
described in the lab frame by, 

z c = 



Vc 



A c cos(wt) 
A c sin(urf) 



(1) 



Also in the lab frame, the components of the vector n 
are given by, 



n x = sin(7) cos(cj£ — a xz ) 
n y = sin(7) sin(ut - a xz ) 
n z = cos(7) 



(2) 



where a xz is the phase shift between the circular motion 
of the pan center, and the precessional motion of n. 

In order to analyze the behavior of the granular parti- 
cles in the moving pan, we also need to derive expressions 
describing the motion, in the lab frame, of an arbitrary 
point P fixed on the pan bottom. Let {x p ,y p ) be the 
Cartesian coordinates of P in a frame of reference fixed 
to, and lying in the plane of, the pan bottom. The co- 
ordinates of P in the lab frame, (x,y,z), are then given 
by, 



x(t) 

y(t) 

z(t) 



x p + A c cos(ujt) 
y p + A c sin(wt) 



(3) 
(4) 

(5) 



The expression for z{t) above is determined from the 
condition that the pan is a rigid body, and from the 
orthogonality of the vectors n and p, where p is the 
vector from the center of the pan bottom to the point 
P, in the pan frame. If we express (x p ,y p ) in terms 
of a polar coordinate system (r, 9) that is also fixed 
in the pan frame, then x p — r cos 9 and y p = r sin 9, 

where r = \Jx p + y p . Using the trigonometric identity 

cos a cos (3 + sin a sin /3 = cos(a — (3), Eq.[5]can be rewrit- 
ten as, 



z(t) = —r tan(7) cos(wi 



(6) 



Note that in Ref. [ij , A c is referred to as x max = y m ax 
and the maximum vertical amplitude of vibration i?sin7, 
where R is the radius of the pan, is referred to as z max . 
With these changes of notation, the representation of the 
pan motion presented in Eqs.|310]and[f)]is consistent with 
Eqs. 1-3 of reference 



III. NUMERICAL SIMULATIONS 

A. Event Driven Algorithm 

The event driven algorithm [T3 | used in our simulations 
has proved very effective in simulating agitated granular 
matter. For example, simulations of standing wave pat- 
terns in 3D vertically oscillated granular layers [13 show 
remarkable agreement with experiments. The algorithm 
itself is well described elsewhere (see e.g. Refs. |l4i llq L 
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In short, it involves five steps: (i) the calculation of colli- 
sions times between particles and between particles and 
walls; (ii) determining which pair of objects are next to 
collide by finding the minimum of the collision times; (iii) 
evolving the system forward in time to the next-to-occur 
collision by changing the coordinates and velocities of all 
particles according to the equations of motion of free par- 
ticles in a gravitational field; (iv) changing the velocities 
of the two colliding particles (or of a particle hitting a 
wall) according to the conservation of linear and angular 
momentum and a specified rule for the dissipation of en- 
ergy in the collision; and (v), returning to step (i). The 
algorithm assumes that interactions occur only through 
instantaneous binary collisions. 

We have adapted this algorithm to our specific system 
by implementing the interactions between particles and 
moving walls of a cylindrical flat-bottomed container in a 
manner that is consistent with hard sphere dynamics. As 
described below, we give special care to the implementa- 
tion of roughness of the bottom surface of the container. 
We have tested several models for the interaction between 
particles and a rough bottom plate and decide on the one 
presented below because of its simplicity combined with 
its efficiency. 



B. System Parameters 

In all our simulations the particles are monodisperse 
spheres of diameter a possessing both linear and angular 
momentum. Otherwise, our numerical model has a fairly 
large number of parameters that have to be set so that the 
model can be validated against the experimental results 
of Ref . 1] . In the following we organize these into three 
groups: dissipative, system size, and forcing parameters. 

Dissipative Parameters: These parameters control the 
way energy is dissipated in the system. We implement 
inelastic collisions between particles in the manner de- 
scribed in . The collisions between particles are thus 
dependent on the coefficient of restitution e, defined as 
the ratio of the normal components of the relative veloc- 
ities of two colliding particles before and after the col- 
lision; the rotational coefficient of restitution (3, defined 
as the ratio of tangential components of relative surface 
velocities of colliding particles before and after the col- 
lision; and the coefficient of friction /1 for the material 
from which the particles are made. 

We set e to be velocity dependent in the following way: 

e(«n) = ( 1 - Bvl ' A ' Vn < v ° (7) 

Here v n is the component of relative velocity along the 
line joining particle centers, B = (1 — e)v , v a = ^fgd 
is a characteristic velocity of the particles, g is the accel- 
eration due to gravity, and e is a tunable parameter of the 
simulations. Decreasing e increases the energy dissipated 
during each collision. The advantages of this model for 



the coefficient of restitution are discussed and demon- 
strated in Refs. El El El El EH El. We have set 
e = 0.6, because this value gives the best agreement in 
terms of particle compaction observed in similarly agi- 
tated states in experiments with nylon balls and in our 
simulations. 

The dissipation of rotational kinetic energy due to 
colliding particles is controlled by (3. As validated in 
Ref. , the dependence of (3 on the relative particle ve- 
locities and on /1 depends on the value of the threshold 
parameter j3 Q that controls the choice between a rolling 
or sliding solution in a given particle-particle collisions. 
Ref. E3 suggests that (3 a = 0.35 is a reasonable choice for 
plastic sphere collisions based on comparison with exper- 
iments. We use this value in all our simulations. We have 
set (i = 0.5 in all our simulations based on the findings 
of Ref. 

The dissipation of energy in collisions between parti- 
cles and the side walls and the bottom of the container is 
determined by the same set of parameters (e, /i, and (3 ) 
as in particle-particle collisions. For the case of collisions 
between particles and the side walls we have observed, 
however, that the simulated system's behavior is not de- 
tectably changed when using e — 1 (elastic collisions) 
instead of e < 1 (inelastic collisions). Thus, in all our 
simulations we set e = 1 for collisions between particles 
and the side walls. 

Special attention is given to the implementation of 
momentum transfer and dissipation of energy in colli- 
sions between particles and the rough bottom of the 
container. Instead of using three dissipative parame- 
ters as in particle-particle collisions, we implement a one- 
parameter model as follows. In the reference frame of the 
moving pan, the particle strikes the bottom plate with a 
relative incident velocity u at the contact point where 
the surface normal is n. The particle is reflected from 
the surface with relative velocity v — e ( | u | ) [u — ( 2u ■ n) n] 
where e(|u|) has the same form as Eq. except that it 
is a function of the absolute value of the incident veloc- 
ity and not just the normal component. Thus, if not for 
the factor of e(|u|), the particle-bottom collision would 
be elastic. The decrease of all three components of the 
reflected velocity by the factor e(|u|) is essential here for 
modeling the entrainment of particles by the rough bot- 
tom surface. Indeed, the effect of this decrease is that 
the particle's velocity (in the lab frame) after the colli- 
sion tends to be more similar to the velocity of the rough 
bottom surface at the point of contact, than before the 
collision. This, in turn, makes the choice of the reflection 
direction of little importance. Statistically the average 
direction of reflection must be in the plane of the vectors 
u and n, as it is in our model. Thus the only parameter 
of our model for the rough bottom surface is e in Eq. [7| 
We have tested several values of this parameter in the 
range between 0.6 and 0.9. Lower values of e provide 
better entrainment of particles by the bottom surface, 
but at the same time make the average collision time 
close to the surface so small that roundoff errors ren- 
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der the simulations unphysical. Thus, we have adopted a 
value of e = 0.8 as a compromise that realizes particle en- 
trainment at the bottom surface as well as well-behaved 
numerical simulations. 



System Size Parameters: These parameters include the 
particle diameter a, number of particles N and the radius 
of pan R. The layer depth H can be estimated assum- 
ing the volume fraction $ = (N^ira 3 )/(irR 2 H) is equal 
to the experimentally determined 15] value of 0.58. For 
example, in the laboratory experiments of Ref. Q using 
nylon balls, a 30 kg load of particles (each weighing 0.9 g 
and 11.5 mm in diameter) in a pan of radius 0.381 m 
gave H/a = 8.7 for the layer depth. In the present simu- 
lations we set the number of particles to either N = 2175 
or N = 15000. In the first case we have R/a = 12.5 and 
H/a = 4.0, and in the second case we have R/a = 33 
and also H/a — 4.0. Although the layer depth in our 
simulations is smaller than in the experiments, as we will 
show below, we are able to reproduce the main features of 
the dynamical behavior observed in experiments. In this 
work we have chosen to present the results for only the 
smaller system size (N — 2175) because the shorter com- 
puting times allow for larger statistical sampling. Our 
tests of the system behavior using the larger number of 
particles indicates that the results are qualitatively the 
same as for the smaller system. 



Forcing Parameters: The frequency of vibrations / = 
u/2tt, the phase shift a xz , the amplitude of horizontal 
vibrations A c , and the tilt angle 7, were introduced in 
Section [H] For all results presented here we consider the 
set (/, a xz ) as the independent variables and we study the 
system behavior as a function of this set, with all other 
parameters fixed. Table |U shows the values (or range of 
values) of the forcing and system size parameters studied 
here, compared to those examined in the experiments of 
Ref. 0. 



As stated above, we conduct simulations on a system 
that is smaller than that of the experiments in order to 
allow us to study as wide a range of dynamical states as 
possible. Since the pan radius R in our simulations is 
therefore smaller than in the experiments, we have set 
the value of the tilt angle 7 to a larger value, such that 
the maximum amplitude of vertical motion at the edge 
of the pan z max = R tan 7 is approximately the same 
as in the experiments. We also find that we must set 
the value of A c , the amplitude of the horizontal motion, 
to a larger value (compared to experiments) in order to 
achieve a robust entrainment of the particles near the 
rough bottom surface. This suggests that the model of 
surface roughness used in our simulations is less effective 
in entraining particle motion than the real system. 



TABLE I: Comparison of simulation parameters used in the 
present work, and experimental parameters from Ref. 0. 





/ (Hz) 




A c /a 


7 


R/a 


H/a 


exp 


10 - 20 


0° - 150° 


0.09 - 0.17 


0.15° - 0.60° 


33.1 


8.7 


sim 


5 - 55 


0° - 360° 


0.5 


1.5° 


12.5 


4.0 
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FIG. 2: Velocity field for toroidal motion at / = 25.5 Hz 
and a xz — 144°. Top panel: the velocities of particles in the 
layer nearest the bottom surface, viewed from above. Each 
arrow length is proportional to the magnitude of the average 
velocity corresponding to each grid location. Note that the 
pan center always rotates in the counter-clockwise direction 
in our simulations, but that under the conditions of the state 
shown here, the tangential flow is clockwise. Bottom panel: 
the velocity field of particles of the same system in a vertical 
cross section passing through the pan center. Here, all ar- 
row lengths are fixed to a constant, and do not represent the 
magnitude of the average velocities. r and h approximately 
locate the center of circulation of the radial flow. 



IV. 



RESULTS 



A. Dynamical modes 

Using the set of parameters described in the previous 
section, we study the behavior of the system as a function 
of / and a xz . We conduct simulations for 15 different fre- 
quencies, in the range between 5 and 55 Hz. For each / 
we set the phase angle a xz to at least 20 different values 
in the range — 360°. We thereby analyze a total of more 
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then 300 state points. In the range of / considered here, 
the motion settles into a steady state in less than 1.0 s of 
physical simulation time, which corresponds to approxi- 
mately 10-50 revolutions of the pan, depending on /. To 
ensure equilibration, we therefore run each simulation for 
at least 2 s before starting to accumulate averages. 

From a visual examination of animations of the steady- 
state behavior, we find that at fixed a xz , three main 
regimes of dynamical behavior occur as / increases: (i) 
At small /, we observe a solid-like heap state that un- 
dergoes a collective tangential rotation about the vertical 
axis. This is qualitatively the same behavior observed at 
small / in Ref. [Jj. (ii) At larger /, typically in the range 
of 15-20 Hz, a transition occurs to a "toroidal" motion, 
having simultaneous tangential and radial circulations, 
which is also very similar to that observed experimen- 
tally, (iii) At large /, the toroidal state crosses over to 
a gas-like state in which the particles have a sufficiently 
large kinetic energy that gravity cannot maintain the par- 
ticles in a dense bed at the bottom of the pan; in this 
regime, the particles spread out over the entire system 
volume. This crossover occurs in the range of / between 
30 and 48 Hz, depending on a xz . 

The simultaneous tangential and radial circulations oc- 
curring in the toroidal regime are shown in Fig. The 
top panel of Fig. [21 shows the velocity field for the par- 
ticles in a horizontal layer that is adjacent to the pan 
bottom. The bottom panel shows the velocity field for 
particles in a vertical plane passing through the center of 
the pan. In these vector field plots we display the velocity 
vectors within the corresponding volume of space aver- 
aged over a local volume of c 3 and over 400 snapshots 
of the system taken 0.001 s apart. Therefore, the total 
averaging time (0.4 s) is much longer than the period of 
oscillation of the pan in the range of / studied here. 



B. Phase diagram 

We next seek to create a dynamical phase diagram of 
the system, i.e. to subdivide the plane of / and a xz 
into domains according to where each of the observed 
dynamical modes occurs. To do so, we quantify the state 
of circulation in the granular bed using two measures of 
the average velocity of the particles. 

The first measure, Ut, quantifies the tangential circula- 
tion and is the average angular velocity of particles in the 
horizontal plane around the vertical line of symmetry: 



k x ?i) 



(8) 



Here, Vti is the tangential component of Vj, the velocity 
of particle i, about the vertical axis, di is the perpen- 
dicular distance from a particle i to the z-axis of the lab 
frame; k is the unit vector pointing along the z-axis of the 
lab frame; and f , is the radial unit vector at the position 
of particle i in the lab frame. 
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FIG. 3: uj r (triangles) and uo t (circles) as a function of a xz at 
various /. To facilitate comparison, uj r values are multiplied 
by 30, and u t values by 1.5. 



The second measure quantifies the radial circulation. 
We define this measure, u> r , as the average angular veloc- 
ity of particles about the curved axis that forms a hori- 
zontal circle of radius r D and elevated above the bottom of 
the pan at a distance h a . As illustrated in the lower panel 
of Fig. [21 if we choose r = 0.75R and h a — 0.75H, we 
find that this curved axis approximately coincides with 
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FIG. 4: Phase diagram of dynamical modes found in sim- 
ulations. The solid-like heap state occurs at low /, in the 
diagonally striped areas. The domain of the gas-like state at 
high / is denoted by the square lattice pattern. The remain- 
ing white areas, bounded by triangles, indicate the domains in 
which the toroidal circulation is observed. The label "inward" 
indicates the region where the radial flow is such that particles 
travel toward the center along the top surface of the bed. The 
label "outward" indicates the region where the radial flow is 
such that particles travel away from the center along the top 
surface of the bed. For both the heap and toroidal states, 
the line that connects the open circles indicates a change of 
direction of the tangential flow. Left (right) arrows indicate 
that the direction of the tangential flow is against (with) the 
direction of the pan rotation. The dashed rectangle shows 
the approximate region corresponding to the experiments of 
Ref. 0]. Error bars for the boundaries of the toroidal domains 
(triangles) are of the size of the symbol for low / and increase 
to up to three times the size of the symbol for high /. 



V 1 

/_ . \ \ 

^ ~- N \ \ \ 

„ _ \ \ \ I I I 
/III/// 



\ I 



I / / / / ^ ^ 

!!////_ 
\ \ \ \ \ I ^ 



s 1 / ^ / ^\ \ \ \ \ 



FIG. 5: Velocity field in a vertical cross section passing 
through the pan center, illustrating outward radial flow, oc- 
curring at / = 27 Hz and a xz = 0° . 



the center of the radial circulation observed in our simu- 
lations. uj r is then given by, 



1 V r ^i 



1 

N 2^' 



1 



Vj • (tj x aj) 



(9) 



Here i> r ,i is the tangential component of v i; in a vertical 
plane, with respect the circulation center defined by r Q 
and h ] and t, = (k x fj). and a, are respectively the 
direction and magnitude of the vector a^ = — (/i D k + 



r r i) which defines the position of particle i with respect 
to the center of circulation in the vertical plane. 

Examples of the behavior of oj t and u> r , as a function of 
a xz at various fixed /, are shown in Fig. 01 I n Fig-GIa), 
we see that cu r = for a xz < 90°, becomes non-zero for 
higher a xz , but then returns to zero as a xz approaches 
360°. Based on our visual observations of the computer 
animations, the heap state exhibits no radial circulation, 
and so the progression of behavior of oj r in Fig. 0{a) re- 
flects the transition first from the heap to the toroidal 
mode, and then back again, as a function of a xz . On the 
other hand, we see that Lu t varies smoothly as a function 
of a xz in Fig.[3Ia), reflecting the fact that tangential cir- 
culation is common to both the heap and toroidal modes. 
Based on these observations, for the purpose of construct- 
ing a phase diagram, we define the toroidal state as that 
in which u> r is clearly different from zero. We also note 
that the variation of tut and u> r with a xz is approximately 
sinusoidal, with the behavior of to r lagging that of uj t by 
about 90°. We will return to this behavior in the next 
section. 

The resulting phase diagram is shown in Fig. 0] The 
most notable behavior of the observations summarized 
in Figs. |21 and 0| is that we find that both the radial and 
tangential circulations reverse their direction as a func- 
tion (primarily) of a xz , as indicated by the occurrence of 
both positive and negative values of u) t and u> r in Fig. [3] 
The result is a rich and highly controllable spectrum of 
dynamical behavior: at low /, we see heap states that ro- 
tate either clockwise or counter-clockwise (viewed from 
above); at higher /, toroidal states with any combination 
of clockwise or counter-clockwise tangential circulation, 
and inward or outward radial circulation (based on the 
flow direction on the upper free surface) can be realized. 
The velocity field of a system with an outward radial cir- 
culation is shown in Fig. [S] In addition, the magnitudes 
of both tut and uj r can be readily tuned through their 
dependence on / and a xz . 

The phase behavior shown in Fig. 0] is in quantitative 
agreement with the experimental results of Ref. 0] for 
the boundary between the heap and toriodal modes. The 
forcing conditions studied in Ref. 1] correspond approx- 
imately to 10 < / < 20 and 100° < a xz < 190°, in the 
notation of this work. In this range, the tangential circu- 
lation is (almost always) in the same direction for both 
the heap and toroidal modes, and, as in experiments, is 
opposite to the direction the pan center rotates. The ex- 
ception, in both simulations and experiments, occurs for 
small values of both / and a xz , where the tangential cir- 
culation occurs in the same direction as the pan center 
rotates. Since the experiments only studied a small frac- 
tion of the phase diagram explored here, the reversal of 
the radial circulation was not observed in experiments, 
but is revealed here, along with the complete demarcation 
of the reversal boundary for the tangential circulation. 

We note that the more subtle modifications of the 
toroidal mode observed experimentally (the "surface 
waves" and "sectors" described in Ref. [l|) are not ob- 
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served here, suggesting that these phenomena depend on 
aspects of the motion and/or particle interactions that 
are not modeled in the present simulations. One pos- 
sibility is that these modifications of the toroidal mode 
are due to variations of the pan motion from a perfect 
sinusoid. As indicated in Ref. the experimental appa- 
ratus was not perfectly symmetrical, and the measured 
motion indicated the presence of a weak pattern of beats 
superimposed on the principal motion. While it is possi- 
ble that the experimentally observed surface waves and 
sectors had there origin in these additional components 
of the motion, more research, beyond the scope of that 
described here, would be required to clarify this issue. 

V. ORIGINS OF OBSERVED FLOWS 

In this section, we discuss some of the physical mech- 
anisms that may underlie the behavior noted so far. We 
focus in turn on two questions: (i) Why do the horizontal 
and radial circulations reverse their direction as a func- 
tion of a xz l (ii) What parameters determine the onset 
of heap-to-toroidal flow? 

A. Circulation reversals 

In order to clarify the nature of the interaction between 
the granular bed and the moving pan, we first examine 
the nature of the contact between the two as a function of 
time. To determine where and when the bed is in contact 
with the pan as a function of time, we subdivide the pan 
bottom into a grid of squares of area a 2 . Over a time 
interval of 0.05 of the period of revolution of the pan, 
we then find the average change of momentum of parti- 
cles colliding with the pan bottom in each grid square. 
This provides a measure of the instantaneous pressure 
exerted by the bed on each area element of the pan bot- 
tom. When the system is in a dynamical steady state, 
we can average over "snapshots" of the pressure distri- 
bution on the pan bottom at different times by rotating 
each snapshot around the pan center such that the di- 
rection given by 9 = tot always points along the positive 
x-axis. Pressure contours resulting from this averaging 
process over a total time of 300 revolutions of the pan 
are shown in Fig.[S]for the case / = 23.4 Hz, a xz = 144°. 
The pressure contours show that a well-defined "zone of 
contact" exists between the granular bed and the pan 
bottom at a particular angular position. Similar results 
are obtained at other values of / and a xz , and differ only 
in the angular position of the zone of contact. 

We find that 9 , the angular position of the center of 
the zone of contact is, in all cases examined, close to the 
value, 

6 = ujt- a xz - tt/2 (10) 

A relatively simple rationale supports this expression: 
We expect that the contact between the bed and the 




FIG. 6: Interrelationship of the zone of contact between the 
granular bed and the pan bottom, and the motion of the pan, 
as viewed from above, for / = 23.4 Hz, a xz = 144°. The con- 
tours are lines along which the average momentum change of 
particles colliding with the bottom surface is constant. The 
vector labeled n is the projection of the pan normal vector 
into the horizontal plane, and indicates the direction along 
which the vertical displacement of the pan is most negative. 
The vector labelled — n therefore indicates the direction along 
which the vertical displacement of the pan is most positive. 
Since the pan center orbits counter-clockwise in the horizon- 
tal plane, the half circle above and to the left of the bisector 
formed by n and — n is that half of the bottom surface that 
is rising. The zone of contact, located in the region around 
the center of the momentum-change contours, is located on 
the rising side of the pan bottom, and has its center near the 
direction o = cot — a xz — n/2. The horizontal "pushing veloc- 
ity" of the bottom pan surface, v p , and its polar component 
vectors, are also shown. 

pan occurs on that side of the pan which is rising from 
its maximum negative displacement in the z direction 
to its maximum positive displacement. As shown in 
Fig. El when the normal vector of the pan n is precessing 
counter-clockwise, the "rising" half of the pan bottom is 
identified by those values of 9 such that, 

< (wt - a xz -9)<tt (11) 

We see from Fig. \G\ that the maximum pressure of the 
zone of contact is approximately in the middle of this 
half of the pan bottom, the angular position of which is 
given approximately by Eq. 1101 

Next we note that we have modeled the pan bottom 
as a rough surface that entrains the motion of particles 
that collide with it; i.e. particles tend to have a post- 
collision velocity with a direction similar to that of the 
surface element of the pan bottom that they strike. Since 
the contact between the granular bed and the pan bot- 
tom is localized to the zone of contact, we expect that 
the velocity of the surface elements at the angular po- 
sition given by Eq. 1101 will have a significant influence 
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on the direction of flow of the bed. Let us define the 
"pushing velocity" v p as the component of the pan bot- 
tom velocity in the horizontal plane at 9 = o . The time 
derivatives of Eqs. and 0] give us the x and y compo- 
nents of v p : Vp tX = — A c lu sin((jjt) and v PjV — A c to cos(cut). 
Converting these Cartesian components to those in a 
polar coordinate system gives the radial and tangen- 
tial components of v p . The desired radial unit vector 
is f = (cos6> , sm0 o ) and the tangential unit vector is 
t = [cos(0 o + 7r/2),sin(0 o + 7r/2)]. The radial and tan- 
gential components of v p are therefore given by, 

Vp. r = Vp-f = -A c u> cos(a xz ) (12) 
v p ,t = Vp-t = -A c uj sin(a xz ) (13) 

Comparing these expressions with the dependence of 
u) t a- n d w r on a xz shown in Fig. [31 shows that there is a 
good correspondence between the polar components of v p 
and the observed angular velocity of both the radial and 
tangential flows. uj t shows an especially good match to a 
function proportional to — sh^a^z). The correspondence 
of Lu r to — cos(a xz ) is complicated by the fact that uj r = 
in the heap state. However the observed trend for uo r to 
lag u>t by 90° is consistent with the implications of Eas.lT2l 
and [T2I 

Overall, we therefore find that the assumption that the 
particles are entrained by the pan bottom in a specific 
zone of contact explains well the observed circulation. If 
the polar components of v p drive the motion observed 
via u>t and u> r , then Eas. 1121 and 1131 predict the occur- 
rence of both clockwise and counter-clockwise tangential 
circulations, as well as both inward and outward radial 
flows. In particular, Eas. 1121 and 1131 also provide a basis 
for understanding why the tangential flow can be both 
the same as, or opposite to, the direction of rotation of 
the pan center. Most importantly, the dominant role of 
the phase shift a xz in controlling the direction of the ob- 
served flows is confirmed. 

B. Onset of toroidal flow 

Since both the heap and toroidal circulations exhibit 
tangential flow, understanding the onset of the toroidal 
mode requires us to understand the conditions required 
for the appearance of the radial circulation. This is a 
complex question, since the heap-to-toroidal boundary 
in Fig. 0] is clearly dependent both on a xz and /. 

To make a preliminary attempt to understand the de- 
pendence of the heap-to-toroidal boundary on a xz and 
/, let us consider a ri the maximum radial acceleration 
experienced by the particle bed, which drives the radial 
circulation observed in the toroidal mode. By analogy 
with many other systems exhibiting granular circulation, 
we expect that the onset of the toroidal mode occurs 
when a r exceeds a critical threshold large enough to sus- 
tain the radial circulation. 

To estimate the scaling of a r with a xz and /, we first 
assume that a r = A p tu 2 , where A p is the average displace- 
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FIG. 7: Model phase diagram based on Eqs. 1121 1131 and 
1151 Solid lines indicate the onset of radial flow. Dotted lines 
indicate the change of direction of the tangential flow. Arrows 
and region labels have the same meaning as in Fig. 2] 

ment of the particle bed in the radial direction as it is 
pushed by the rough pan bottom. Given the correspon- 
dence between the components of the pushing velocity 
v p and the circulation speeds Wt and ui r established in 
the previous section, we propose that A p is proportional 
to the product of \v v>r \ and t p , the time interval during 
which particles are exposed to the zone of contact. We 
also assume that t p is simply some fraction of the period 
of the pan motion, and so is proportional to 2tt/uj. Com- 
bining these assumptions with the expression for v p ^ r in 
Eq. El we therefore can write, 

a r cx lu 2 | cos a xz | (14) 

If the onset of the heap-to-toroidal transition occurs at 
a critical value of a r = a r ^ c then the corresponding critical 
angular frequency lo c for the onset of radial circulation 
(and therefore toroidal flow) is given by, 

^coc,/ ° r ' c (15) 
V Icosa^zl 

Fig.[7|shows a plot of the critical frequency f c — u) c /2it in 
the f-a xz plane, based on the proportionality expressed 
in Eq. 1151 with the constant of proportionality chosen so 
that the mimimum onset frequency is 15 Hz, as found 
in our simulations. Comparing Figs. 0] and [7] shows that 
the model expressed in Ea. ll5l captures several qualitative 
features of the observed variation of the heap-to-toroidal 
transition: e.g. the occurrence of maxima in the transi- 
tion boundary as a function of a xz ; and the separation 
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FIG. 8: Three types of boundary condition occurring in this 
system: (a) the periodic boundary condition involving the 
tangential flow; (b) the vertical wall obstacle encountered by 
the inward radial flow; (c) central "plume" obstacle encoun- 
tered by the outward radial flow. 

of the toroidal regime into two distinct domains, one ex- 
hibiting an inward radial flow, and one outward. 

Clearly the approximations of this simple model should 
be refined if it is to have any quantitative predictive 
value. However, it does support the view that entrain- 
ment of the particle bed by the rough bottom is a key 
factor in determining the overall behavior of the system. 

VI. DISCUSSION AND CONCLUSION 

In summary, we have established the following picture 
for the dynamical behavior observed in our simulations, 
and to a large degree, in the experiments of Ref. 0: 
The granular bed is excited by a rough horizontal sur- 
face that is undergoing both vertical and horizontal vi- 
brations. The vertical vibrations impose a periodicity on 
the contact between the bed and the surface, and the 
horizontal vibrations push the granular bed only dur- 
ing the interval of contact. The phase shift relating the 
vertical and the horizontal components of the vibration 
determines the direction and strength of the horizontal 
pushing force applied by the rough surface on the bed. 
The result is a highly tunable spectrum of behavior, in 
which both tangential and radial circulations, in arbi- 
trary combinations, can be generated. 

The cylindrical geometry studied here plays an impor- 
tant underlying role in the granular flows we observe. 
The tangential circulation is in many ways similar to 
the kind of linear flow observed in vibratory convey- 



ors [23L |24|. but wrapped around into a circle, creating 
in effect a periodic boundary condition. Indeed, a recent 
work studied the flow of granular matter in an annular 
vibratory conveyor and, as in our work, noted both the 
occurrence of a reversing tangential flow, and the con- 
trolling influence on this flow of the phase shift between 
vertical and horizontal vibrations [Hj • The fact that the 
granular bed in our system encounters no obstacles as it 
moves about the center accounts for the ease with which 
the tangential flow occurs not only in the toroidal mode, 
but also the heap mode at very small /. 

In contrast, the granular flow in the radial direction en- 
counters a vertical barrier that it must surmount: either 
the outer wall in the case of the inward radial flow, or the 
central "plume" that is created when the granular parti- 
cles are pushed toward the center in the outward radial 
flow (see Fig. |SJ). Here, the mechanism that drives the 
flow along the bottom surface is also similar to that of a 
vibratory conveyor, but with the complication that this 
driving force must be large enough to overcome the grav- 
itational barrier imposed at the vertical boundary. This 
largely explains why the radial flow only occurs above a 
critical value of /. 

Our results also demonstrate that complex dynamical 
modes observed in real (i.e. industrial) granular matter 
devices can, on the one hand, be understood using com- 
puter simulations; and on the other hand, the behavior 
observed in these real devices can reveal important fun- 
damental principles (such as the role of the phase factor 
established here) relating to the creation and control of 
highly specialized granular flows. 
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